CIND820 Project Course Code

Ryan Boyd

Libraries imported

import laspy as lp, sklearn as skl, numpy as np, matplotlib as mp, pandas as pd

from sklearn import cluster

from sklearn import preprocessing as prep

import matplotlib.pyplot as plt

import PIL
from PIL import ImageStat as istat
from PIL import ImageOps

Data loaded and prepared

path_to_data = "F:/Data/Lidar/dtvan/dtvan.laz"
with lp.open(path_to_data) as las_file:
    las_data = las_file.read()

#print("Total points:" + str(las_data.header.point_count))

Total points:20975490

#print("Classes: " + + str(set(list(las_data.classification))))

{1, 2}

df = pd.DataFrame({"X":list(las_data.x),"Y":list(las_data.y),"Z":list(las_data.z),"Intensity":las_data.intensity,"return_num":las_data.return_number,"totalreturns":las_data.num_returns,"classification":las_data.classification})

fig,ax = plt.subplots(figsize = (15,15))
ax.scatter(df['X'],df['Y'],zorder=1,alpha=0.25,c='black',s=0.001)

<matplotlib.collections.PathCollection at 0x2177a8c99a0>

The full original dataset is too large to work with so it is clipped to a smaller study area which was arbitratily selected. The dataframe is queried by the study area longitude and latitude boundary. The pre-classified ground points are removed (since we are not concerned with ground right now), which is class 2, class 1 is the unclassified points and we only want to work with these.

# Define the area of interest, these values in meters are in the "NAD83 UTM 10" coordinate system of the provided dataset
# These are the upper and lower limits in meters to be used, these can be found using google maps or other free sites/software
# These were selected somewhat arbitrarily -rb
aoi_extent = {'xmax':492349.0731766,'xmin':492043.6935073,'ymax':5458645.8660691,'ymin':5458340.4864470}

#query the dataframe to return only the points within the extent above and remove the points defined as ground as well -rb
df_clip = df.query("X>{0}&X<{1}&Y>{2}&Y<{3}&classification==1".format(aoi_extent['xmin'],aoi_extent['xmax'],aoi_extent['ymin'],aoi_extent['ymax']))

df_clip

	X	Y	Z	Intensity	return_num	totalreturns	classification
4394229	492044.02	5458475.09	18.75	27	2	2	1
4394231	492044.00	5458475.08	18.74	26	2	2	1
4394233	492043.98	5458475.08	18.75	26	2	2	1
4394235	492043.94	5458475.08	18.78	32	2	2	1
4394237	492043.91	5458475.08	18.75	41	2	2	1
...	...	...	...	...	...	...	...
20825059	492345.70	5458600.09	15.56	153	1	1	1
20825060	492345.82	5458600.05	14.90	116	1	1	1
20825061	492346.11	5458600.05	14.84	139	1	1	1
20825062	492346.34	5458600.04	14.50	66	1	2	1
20825063	492346.61	5458600.04	14.41	29	1	2	1

789354 rows × 7 columns

# Display the points on a scatter plot -rb
ig,ax = plt.subplots(figsize = (15,15))
ax.scatter(df_clip['X'],df_clip['Y'],zorder=1,alpha=0.5,c='black',s=0.01)

<matplotlib.collections.PathCollection at 0x21757766490>

# And we can use the height value to visualize in 3d, since the values are in m for all 3 axes, it plots nicely as is -rb
fig,ax = plt.subplots(figsize = (15,15)),plt.axes(projection='3d')
ax.scatter3D(df_clip['X'],df_clip['Y'],df_clip['Z'],c='black',s=0.01,alpha=0.5)
#from matplotlib import cm
#ax.plot_surface(df_clip['X'],df_clip['Y'],df_clip['Z'],cmap=cm.coolwarm,linewidth=0,antialiased=False)

<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x2178b6284f0>

df_ground = df.query("X>{0}&X<{1}&Y>{2}&Y<{3}&classification==2".format(aoi_extent['xmin'],aoi_extent['xmax'],aoi_extent['ymin'],aoi_extent['ymax']))

df_ground['Z']
# The Z values are an absolute elevation which is not as useful as height relative to the ground
# I need to create a ground image to subtract from the elevation in order to get the height of the object relative to the ground -rb

6437554     5.20
6441448     5.39
6441451     5.25
6441453     5.24
6441456     5.26
            ... 
20822410    5.51
20822411    5.49
20822412    5.50
20822414    5.52
20822415    5.53
Name: Z, Length: 75058, dtype: float64

# Plotting the points that were classified as ground in the provided dataset
fig,ax = plt.subplots(figsize = (10,10))
ax.scatter(df_ground['X'],df_ground['Y'],c='black',s=0.01,alpha=0.5)

<matplotlib.collections.PathCollection at 0x21778b40a00>

Data normalized and preprocessed for analysis

x_normal = (df_clip['X'] - min(df_clip['X']))/(max(df_clip['X']-min(df_clip['X'])))

y_normal = (df_clip['Y'] - min(df_clip['Y']))/(max(df_clip['Y']-min(df_clip['Y'])))

z_normal = (df_clip['Z'] - min(df_clip['Z']))/(max(df_clip['Z']-min(df_clip['Z'])))

i_normal = (df_clip['Intensity'] - min(df_clip['Intensity']))/(max(df_clip['Intensity']-min(df_clip['Intensity'])))

df_normal = pd.DataFrame({'X':x_normal,'Y':y_normal,'Z':z_normal,'I':i_normal,'return_num':df_clip['return_num'],'total_returns':df_clip['totalreturns']})

df_normal

	X	Y	Z	I	return_num	total_returns
4394229	0.001048	0.440777	0.241716	0.006407	2	2
4394231	0.000982	0.440744	0.241629	0.006161	2	2
4394233	0.000917	0.440744	0.241716	0.006161	2	2
4394235	0.000786	0.440744	0.241974	0.007639	2	2
4394237	0.000688	0.440744	0.241716	0.009857	2	2
...	...	...	...	...	...	...
20825059	0.988964	0.850116	0.214187	0.037457	1	1
20825060	0.989357	0.849985	0.208492	0.028339	1	1
20825061	0.990307	0.849985	0.207974	0.034007	1	1
20825062	0.991060	0.849953	0.205040	0.016018	1	2
20825063	0.991944	0.849953	0.204263	0.006900	1	2

789354 rows × 6 columns

# Plotting normalized looks the same but with new scale
fig,ax = plt.subplots(figsize = (10,10))
ax.scatter(df_normal['X'],df_normal['Y'],c='black',s=0.01,alpha=0.5)

<matplotlib.collections.PathCollection at 0x21778e7c6d0>

Add Imagery Data

2015 Imagery data was obtained from the City of Vancouver to extract the RGB values

The image was clipped using external software (QGIS, open-source mapping program) to the same area of interest as above

The selected image size is 4084x4084, the lidar data is normalized by 4084 to extract the nearest pixel value(r,g,b) from the image for each point

image_path = "F:/Data/Lidar/images/BCVANC15_P9_aoiclip.tif"

img = PIL.Image.open(image_path)

rgb_img = img.convert("RGB")

rgb_img

#this can be used to crop the imagery if we knew the exact coordinates, but I used QGIS to clip the imagery instead
#left,top,right,bottom = 0,0,4084,4084
#rgb_img = img.crop((left,top,right,bottom))

#import math
#math.sqrt(istat.Stat(rgb_img).count[0])

# This is needed to align the pixels with the lidar points even though it isn't a perfect fit, it still works approximately
img.size

(4084, 4084)

#The image origin (top left) is different than the coordinate system of the lidar so the image needs to be flipped for the calculation to align them
rgb_img = PIL.ImageOps.flip(rgb_img)

#rgb_img.getpixel((2070,2070))

x_img = (df_clip['X'] - min(df_clip['X']))/(max(df_clip['X']-min(df_clip['X'])))*4083
y_img = (df_clip['Y'] - min(df_clip['Y']))/(max(df_clip['Y']-min(df_clip['Y'])))*4083

y_img

4394229     1799.691522
4394231     1799.557815
4394233     1799.557815
4394235     1799.557815
4394237     1799.557815
               ...     
20825059    3471.024659
20825060    3470.489832
20825061    3470.489832
20825062    3470.356125
20825063    3470.356125
Name: Y, Length: 789354, dtype: float64

coord_array = np.array(pd.DataFrame({"X":x_img,"Y":y_img}))

coord_array

array([[4.27861283e+00, 1.79969152e+03],
       [4.01119953e+00, 1.79955782e+03],
       [3.74378623e+00, 1.79955782e+03],
       ...,
       [4.04342283e+03, 3.47048983e+03],
       [4.04649808e+03, 3.47035613e+03],
       [4.05010816e+03, 3.47035613e+03]])

The nearest R,G,B pixel from the image is extracted and the coordinates are saved as a field to the data frame to be used in the model

rgb_data = []
rgb_data_r = []
rgb_data_g = []
rgb_data_b = []
for coord in coord_array:
    rgb=rgb_img.getpixel((coord[0],coord[1]))
    r=rgb[0]
    g=rgb[1]
    b=rgb[2]
    rgb_data.append(rgb)
    rgb_data_r.append(r)
    rgb_data_g.append(g)
    rgb_data_b.append(b)
df_normal['rgb'] = rgb_data
df_normal['r'] = rgb_data_r
df_normal['g'] = rgb_data_g
df_normal['b'] = rgb_data_b

df_normal

	X	Y	Z	I	return_num	total_returns	rgb	r	g	b
4394229	0.001048	0.440777	0.241716	0.006407	2	2	(0, 0, 0)	0	0	0
4394231	0.000982	0.440744	0.241629	0.006161	2	2	(0, 0, 0)	0	0	0
4394233	0.000917	0.440744	0.241716	0.006161	2	2	(0, 0, 0)	0	0	0
4394235	0.000786	0.440744	0.241974	0.007639	2	2	(0, 0, 0)	0	0	0
4394237	0.000688	0.440744	0.241716	0.009857	2	2	(0, 0, 0)	0	0	0
...	...	...	...	...	...	...	...	...	...	...
20825059	0.988964	0.850116	0.214187	0.037457	1	1	(132, 129, 105)	132	129	105
20825060	0.989357	0.849985	0.208492	0.028339	1	1	(136, 131, 105)	136	131	105
20825061	0.990307	0.849985	0.207974	0.034007	1	1	(135, 130, 101)	135	130	101
20825062	0.991060	0.849953	0.205040	0.016018	1	2	(138, 134, 102)	138	134	102
20825063	0.991944	0.849953	0.204263	0.006900	1	2	(142, 135, 107)	142	135	107

789354 rows × 10 columns

The R,G,B values are normalized like the rest:

r_normal = (df_normal['r'] - min(df_normal['r']))/(max(df_normal['r']-min(df_normal['r'])))
g_normal = (df_normal['g'] - min(df_normal['g']))/(max(df_normal['g']-min(df_normal['g'])))
b_normal = (df_normal['b'] - min(df_normal['b']))/(max(df_normal['b']-min(df_normal['b'])))
df_normal['r'] = r_normal
df_normal['g'] = g_normal
df_normal['b'] = b_normal

df_normal

	X	Y	Z	I	return_num	total_returns	rgb	r	g	b
4394229	0.001048	0.440777	0.241716	0.006407	2	2	(0, 0, 0)	0.000000	0.000000	0.000000
4394231	0.000982	0.440744	0.241629	0.006161	2	2	(0, 0, 0)	0.000000	0.000000	0.000000
4394233	0.000917	0.440744	0.241716	0.006161	2	2	(0, 0, 0)	0.000000	0.000000	0.000000
4394235	0.000786	0.440744	0.241974	0.007639	2	2	(0, 0, 0)	0.000000	0.000000	0.000000
4394237	0.000688	0.440744	0.241716	0.009857	2	2	(0, 0, 0)	0.000000	0.000000	0.000000
...	...	...	...	...	...	...	...	...	...	...
20825059	0.988964	0.850116	0.214187	0.037457	1	1	(132, 129, 105)	0.517647	0.505882	0.413386
20825060	0.989357	0.849985	0.208492	0.028339	1	1	(136, 131, 105)	0.533333	0.513725	0.413386
20825061	0.990307	0.849985	0.207974	0.034007	1	1	(135, 130, 101)	0.529412	0.509804	0.397638
20825062	0.991060	0.849953	0.205040	0.016018	1	2	(138, 134, 102)	0.541176	0.525490	0.401575
20825063	0.991944	0.849953	0.204263	0.006900	1	2	(142, 135, 107)	0.556863	0.529412	0.421260

789354 rows × 10 columns

Dataset statistics and information if needed

#renaming the data frame for clarity
data = df_normal

mp.pyplot.hist(data['X'])

(array([87314., 86999., 79577., 81244., 66140., 75254., 78820., 76905.,
        73739., 83362.]),
 array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ]),
 <BarContainer object of 10 artists>)

mp.pyplot.hist(data['I'])

(array([7.84475e+05, 3.08600e+03, 9.94000e+02, 3.39000e+02, 1.41000e+02,
        8.90000e+01, 4.10000e+01, 5.40000e+01, 4.30000e+01, 9.20000e+01]),
 array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ]),
 <BarContainer object of 10 artists>)

Testing and looking for improvements to be made

fig,ax = plt.subplots(figsize = (15,15))
ax.scatter(results['X'],results['Y'],c=results['k_cluster'],s=0.1,alpha=0.5)

Initial Classification (unsupervised) attempt using kmeans clustering

Variables: Height, Intensity, Number of Returns

#features = pd.DataFrame({"Z":df_normal['Z'],"I":df_normal['I'],"R":df_normal['r'],"G":df_normal['g'],"B":df_normal['b']})

features = pd.DataFrame({"R":df_normal['r'],"G":df_normal['g'],"B":df_normal['b'],"Z":df_normal['Z'],"I":df_normal['I']})

features

	R	G	B	Z	I
4394229	0.000000	0.000000	0.000000	0.241716	0.006407
4394231	0.000000	0.000000	0.000000	0.241629	0.006161
4394233	0.000000	0.000000	0.000000	0.241716	0.006161
4394235	0.000000	0.000000	0.000000	0.241974	0.007639
4394237	0.000000	0.000000	0.000000	0.241716	0.009857
...	...	...	...	...	...
20825059	0.517647	0.505882	0.413386	0.214187	0.037457
20825060	0.533333	0.513725	0.413386	0.208492	0.028339
20825061	0.529412	0.509804	0.397638	0.207974	0.034007
20825062	0.541176	0.525490	0.401575	0.205040	0.016018
20825063	0.556863	0.529412	0.421260	0.204263	0.006900

789354 rows × 5 columns

X1 = np.array(features)

X1

array([[0.        , 0.        , 0.        , 0.24171557, 0.0064071 ],
       [0.        , 0.        , 0.        , 0.24162927, 0.00616067],
       [0.        , 0.        , 0.        , 0.24171557, 0.00616067],
       ...,
       [0.52941176, 0.50980392, 0.3976378 , 0.20797377, 0.0340069 ],
       [0.54117647, 0.5254902 , 0.4015748 , 0.2050397 , 0.01601774],
       [0.55686275, 0.52941176, 0.42125984, 0.20426303, 0.00689995]])

kmeancluster = cluster.KMeans(5,init='random',n_init=10)

k_clusters = kmeancluster.fit(X1)

len(k_clusters.labels_)

789354

data['k_cluster'] = k_clusters.labels_

results = data

Visualization

fig,ax = plt.subplots(figsize = (15,15))
ax.scatter(results['X'],results['Y'],c=results['k_cluster'],s=0.1,alpha=0.5)

<matplotlib.collections.PathCollection at 0x2178b5bc5e0>

fig,ax = plt.subplots(figsize = (15,15)),plt.axes(projection='3d')
ax.scatter3D(data['X'],data['Y'],data['Z'],c=data['k_cluster'],s=0.01,alpha=1)

<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x21783fee8e0>

Initial Distance Clustering

df2 = pd.DataFrame({"X":data['X'],"Y":data['Y'],"k_cluster":data['k_cluster']})

X2 = np.array(df2)

X2

array([[1.04790909e-03, 4.40776763e-01, 2.00000000e+00],
       [9.82414775e-04, 4.40744015e-01, 2.00000000e+00],
       [9.16920457e-04, 4.40744015e-01, 2.00000000e+00],
       ...,
       [9.90306841e-01, 8.49985264e-01, 4.00000000e+00],
       [9.91060026e-01, 8.49952517e-01, 3.00000000e+00],
       [9.91944199e-01, 8.49952517e-01, 3.00000000e+00]])

#dbscan = cluster.DBSCAN(eps=2,min_samples=10,algorithm="auto",leaf_size=25,n_jobs=-1)

#db_clusters = dbscan.fit(X2)

#results['db_cluster'] = db_clusters.labels_

#fig,ax = plt.subplots(figsize = (15,15))
#ax.scatter(results['X'],results['Y'],c=results['db_cluster'],s=0.1,alpha=0.5)

#fig,ax = plt.subplots(figsize = (15,15)),plt.axes(projection='3d')
#ax.scatter3D(results['X'],results['Y'],results['Z'],c=results['db_cluster'],s=0.01,alpha=1)

#optics = cluster.OPTICS(min_samples=10,max_eps=2,leaf_size=15)

#op_clusters = optics.fit(X2)

#len(set(db_clusters.labels_))
#results['op_clusters'] = op_clusters.labels_

Analysis of results

Z/elevation has a large effect on the classifier

R,G,B causes the shadows to be classified

Updates